library(rethinking)
library(MASS)
# 10E1
log(0.35/(1-0.35))
## [1] -0.6190392
# 10E2
exp(3.2)/(1+exp(3.2))
## [1] 0.9608343
# 10E3
OR <- exp(1.7)
A log-odds of 1.7 corresponds to odds ratio of 5.4739474. That is 447.5% increase in odds.
When events are counted on different temporal or spatial scales.
In the aggregated form, the likelihood is Binomial, with n being the total number of trials. In the disaggregated form, the likelihood is Bernoulli, with n being 1.
When the predicotr increases by one unit, the order (i.e., log) of the expected value of count number in the Possion distribution increases by 1.7.
Because the linear combination of predictors may extend from negative infinity to positive infinity, while the probability of a certain event happens must be between 0 and 1. The logit link function allows the probability parameter to be modelled as a linear combination of predictors.
In Possion distribution we are modeling the expected value of count, which has to be positive. The log link function ensures that.
It would imply that the mean of a Possion GLM can only be between 0 and 1.
The binomial distribution has maximum entropy when only two events are possible, and the expected number of each event is assumed to be constant. The constrait for the possion distribution is essentially the same.
# trace plot for stan model
plot(m10.4_stan)
# compare precis output
precis(m10.4_map, depth=2)
## Mean StdDev 5.5% 94.5%
## a[1] -0.73 0.27 -1.16 -0.30
## a[2] 6.67 3.61 0.90 12.45
## a[3] -1.03 0.28 -1.48 -0.59
## a[4] -1.03 0.28 -1.48 -0.59
## a[5] -0.73 0.27 -1.16 -0.30
## a[6] 0.21 0.27 -0.21 0.64
## a[7] 1.75 0.38 1.14 2.37
## bp 0.82 0.26 0.40 1.24
## bpC -0.13 0.30 -0.61 0.34
precis(m10.4_stan, depth=2)
## Warning in precis(m10.4_stan, depth = 2): There were 3 divergent iterations during sampling.
## Check the chains (trace plots, n_eff, Rhat) carefully to ensure they are valid.
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a[1] -0.73 0.27 -1.15 -0.27 3692 1
## a[2] 11.03 5.28 3.42 18.34 2075 1
## a[3] -1.05 0.28 -1.49 -0.61 3091 1
## a[4] -1.04 0.29 -1.51 -0.61 3354 1
## a[5] -0.73 0.27 -1.17 -0.31 2449 1
## a[6] 0.23 0.27 -0.19 0.66 3304 1
## a[7] 1.82 0.39 1.20 2.44 3438 1
## bp 0.83 0.27 0.40 1.26 2193 1
## bpC -0.14 0.31 -0.63 0.33 3128 1
# compare marginal posterior distributions
pairs(m10.4_map)
pairs(m10.4_stan)
# fit simpler models
m10.1 <- map(
alist(
pulled_left ~ dbinom( 1 , p ) ,
logit(p) <- a ,
a ~ dnorm(0,10)
) , data=d )
m10.2 <- map(
alist(
pulled_left ~ dbinom( 1 , p ) ,
logit(p) <- a + bp*prosoc_left ,
a ~ dnorm(0,10) ,
bp ~ dnorm(0,10)
) , data=d )
m10.3 <- map(
alist(
pulled_left ~ dbinom( 1 , p ) ,
logit(p) <- a + (bp + bpC*condition)*prosoc_left ,
a ~ dnorm(0,10) ,
bp ~ dnorm(0,10) ,
bpC ~ dnorm(0,10)
) , data=d )
# compare models with WAIC
compare(m10.1, m10.2, m10.3, m10.4_map)
## WAIC pWAIC dWAIC weight SE dSE
## m10.4_map 550.5 15.8 0.0 1 18.60 NA
## m10.2 680.5 2.0 130.1 0 9.41 18.14
## m10.3 682.5 3.1 132.0 0 9.39 18.06
## m10.1 687.9 1.0 137.4 0 7.08 18.92
# trace plot for stan model
plot(m10H3_stan)
# compare precis outputs
precis(m10H3_map)
## Mean StdDev 5.5% 94.5%
## a 0.59 0.66 -0.47 1.65
## bp 4.24 0.90 2.81 5.67
## bv -4.59 0.96 -6.13 -3.06
## ba 1.08 0.53 0.23 1.93
precis(m10H3_stan)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a 0.66 0.69 -0.42 1.77 4046 1
## bp 4.64 1.01 3.06 6.18 3115 1
## bv -5.05 1.07 -6.69 -3.35 2615 1
## ba 1.14 0.54 0.27 1.98 4526 1
pairs(m10H3_map)
pairs(m10H3_stan)
Overall the quadratic approximation is OK, but since the marginal posterior distributions for bp and bv are skewed, the quadratic approximation is not entirely accurate.
# the predicted probability of success and it 89% interval
mu <- link(m10H3_stan)
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu.mean <- apply(mu, 2, mean)
mu.HPDI <- apply(mu, 2, HPDI)
# empty plot frame
plot(0, 0, type="n", xlab="P/A/V", ylab="Probability of Success",
ylim=c(0, 1), xaxt="n", xlim=c(1,8))
axis(1, at=1:8, labels=c("1/1/1", "1/1/0", "1/0/1", "1/0/0",
"0/1/1", "0/1/0", "0/0/1", "0/0/0"))
# plot raw proportion
prop <- d$y/d$n
points(1:8, prop, pch=19)
# plot predicted probability and 89% HPDI
points(1:8-0.1, mu.mean)
points(1:8-0.1, mu.HPDI[1,], pch=3)
points(1:8-0.1, mu.HPDI[2,], pch=3)
# the predicted success count and it 89% interval
y.sim <- sim(m10H3_stan)
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
y.mean <- apply(y.sim, 2, mean)
y.HPDI <- apply(y.sim, 2, HPDI)
# empty plot frame
plot(0, 0, type="n", xlab="P/A/V", ylab="Success Count",
ylim=c(0, 30), xaxt="n", xlim=c(1,8))
axis(1, at=1:8, labels=c("1/1/1", "1/1/0", "1/0/1", "1/0/0",
"0/1/1", "0/1/0", "0/0/1", "0/0/0"))
# plot raw count
points(1:8, d$y, pch=19)
# plot predicted success count and 89% HPDI
points(1:8-0.1, y.mean)
points(1:8-0.1, y.HPDI[1,], pch=3)
points(1:8-0.1, y.HPDI[2,], pch=3)
# check trace plot
plot(m10H3_inter)
compare(m10H3_stan, m10H3_inter)
## WAIC pWAIC dWAIC weight SE dSE
## m10H3_inter 93.8 4.6 0.0 0.92 12.73 NA
## m10H3_stan 98.6 4.0 4.8 0.08 13.34 4.72
precis(m10H3_inter)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a -0.77 1.04 -2.33 0.92 2601 1
## bp 6.58 1.43 4.42 8.87 2232 1
## bv -5.30 1.18 -7.15 -3.48 2894 1
## ba 3.44 1.24 1.51 5.42 2339 1
## bpa -2.97 1.37 -5.15 -0.81 2349 1
# the predicted probability of success and it 89% interval
mu <- link(m10H3_inter)
## [ 100 / 1000 ]
[ 200 / 1000 ]
[ 300 / 1000 ]
[ 400 / 1000 ]
[ 500 / 1000 ]
[ 600 / 1000 ]
[ 700 / 1000 ]
[ 800 / 1000 ]
[ 900 / 1000 ]
[ 1000 / 1000 ]
mu.mean <- apply(mu, 2, mean)
mu.HPDI <- apply(mu, 2, HPDI)
# empty plot frame
plot(0, 0, type="n", xlab="P/A/V", ylab="Probability of Success",
ylim=c(0, 1), xaxt="n", xlim=c(1,8))
axis(1, at=1:8, labels=c("1/1/1", "1/1/0", "1/0/1", "1/0/0",
"0/1/1", "0/1/0", "0/0/1", "0/0/0"))
# plot raw proportion
prop <- d$y/d$n
points(1:8, prop, pch=19)
# plot predicted probability and 89% HPDI
points(1:8-0.1, mu.mean)
points(1:8-0.1, mu.HPDI[1,], pch=3)
points(1:8-0.1, mu.HPDI[2,], pch=3)
Comparing the two models with WAIC shows that the model with the interaction term gets 92% of the Akaike weights, while the model without interaction only gets 8%. So the model with the interaction term is estimated to perform much better in predicting future observations.
# check trace plot
plot(m10H4_stan)
# compare precis output
precis(m10H4_map)
## Mean StdDev 5.5% 94.5%
## a -1.46 0.45 -2.18 -0.74
## bp 0.03 0.01 0.02 0.04
precis(m10H4_stan)
## Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## a -1.54 0.48 -2.24 -0.77 1221 1
## bp 0.03 0.01 0.02 0.04 1227 1
pairs(m10H4_map)
pairs(m10H4_stan)
# plot the expected counts
PCTCOVER.seq <- seq(from=0, to=100, by=2)
pred.data <- data.frame(PCTCOVER=PCTCOVER.seq)
count.sim <- sim(m10H4_stan, data=pred.data, n=8000)
## [ 800 / 8000 ]
[ 1600 / 8000 ]
[ 2400 / 8000 ]
[ 3200 / 8000 ]
[ 4000 / 8000 ]
[ 4800 / 8000 ]
[ 5600 / 8000 ]
[ 6400 / 8000 ]
[ 7200 / 8000 ]
[ 8000 / 8000 ]
count.mean <- apply(count.sim, 2, mean)
count.HPDI <- apply(count.sim, 2, HPDI)
plot(SALAMAN ~ PCTCOVER, data=d, col=rangi2,
xlab="Percentage of ground cover", ylab="Counts of salamanders in each plot", pch=19)
lines(PCTCOVER.seq, count.mean)
shade(count.HPDI, PCTCOVER.seq)
# check trace plot
plot(m10H4.1)
plot(m10H4.2)
plot(m10H4.3)
compare(m10H4.1, m10H4.2, m10H4.3)
## WAIC pWAIC dWAIC weight SE dSE
## m10H4.1 213.5 4.8 0.0 0.87 26.42 NA
## m10H4.2 217.6 7.6 4.1 0.11 27.27 1.44
## m10H4.3 221.9 10.4 8.4 0.01 28.18 2.42
count.ensemble <-ensemble(m10H4.1, m10H4.2, m10H4.3)
count.mean <- apply(count.ensemble$sim, 2, mean)
count.HPDI <- apply(count.ensemble$sim, 2, HPDI)
observed.count <- d$SALAMAN+rnorm(n=nrow(d), mean=0, sd=0.1) # dodge points a bit
plot(count.mean ~ observed.count, xlab="Observed Count", ylab="Predicted Count",
xlim=c(0, 12), ylim=c(0, 12), col=rangi2)
abline(a=0, b=1, lty=2)
arrows(observed.count, count.HPDI[1, ], observed.count, count.HPDI[2, ],
length=0.05, angle=90, code=3)